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We present the first numerical code based on the Galerkin method to integrate the field equations 
of the Bondi problem. The Galerkin method is a spectral method whose main feature is to provide 
high accuracy with moderate computational effort. Several numerical tests were performed to verify 
the issues of convergence, stability and accuracy with promising results. This code opens up several 
possibilities of applications in more general scenarios for studying the evolution of spacetimes with 
gravitational waves. 
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D ■ I. INTRODUCTION 

' In their seminal wor ks Bondi et al[ll, 0| have launched the basis for a detailed description for the dynamics of the 
exterior axisymmetric spacetime of a bounded source undergoing a process of gravitational wave emission. They have 
established the suitable form of the metric, the field equations and also exhibited the set of permissable coordinate 
O ', transformations that preserve the nature of the metric of what is known by the Bondi problem. Their most important 
O^- result was the introduction of the news function that determines the rate of mass loss due to gravitational wave 
^ ] extraction, along with providing an invariant characterization of the presence of gravitational radiation. Nonethe- 
less, another important pioneering aspect of this work was the presentation of the field equations in the scheme of 
' characteristics [il more than ten years before the formalization of the Cauchy approach^ of Einstein equations. 

, Due to the complicated form of the Bondi equations there are no analytical solutions describing the evolution 
^ ■ of a non-stationary axisymmetric spacetime endowed with gravitational waves. Therefore, the main features of the 
C ■ I dynamics of the Bondi problem can only be determined through numerical simulations. The first numerical code 
' applied to integrate the Bondi equations was developed by Issacson, Welling and Winicourfsj, and later improved 
^ by a version that avoided instabilities near the vertexiG]. Other groups have also presented distinct strategies for 
. . constructing their codes as, for instance, using the tetrad formalism^ or combining Cauchy and characteristics 
0^ ' evolution that allowed the extension of the Bondi problem to full axisymmetryQ. A detailed description of the 
, numerical schemes used for the characteristic evolution can be found in Ref. [sj, but the common feature shared by 
all of these codes is that they were constructed using finite difference techniques. 

Spectral methods represent an alternative strategy to solve numerically partial differential equations (PDEs), 
and are based on the idea of approximating the solution as a finite series of a convenient set of analytical basis 
functions. As a consequence, any spectral method transforms a PDE or a system of PDEs into a finite system 
of ordinary differential equations for the coefficients that multiply the basis functions, where it is expected that 
'■ the greater is the number of coefficients, or the number of ordinary differential equations, more accurate is the 
approximate solution. There are two main attractive features of the spectral methods: (i) a considerable economy 
of the computational resources to achieve a given accuracy if compared with the finite difference techniques; (ii) the 
possibility of selecting a coordinate system adapted to the geometry of the problem under consideration allowing an 
exact treatment of pseudo-singularities present in the chosen coordinates. As a result the use of spectral methods has 
been increased considerably in the last years and figures as viable strategy for Numerical Relativity '11]. In particular, 
amongst the several types of spectral methods we are interested in the Galerkin method(l3|. whose main feature is 
that each basis function is selected so as to satisfy automatically the boundary conditions. 

In this work we have developed the first numerical code based on the Galerkin method [13, to integrate the field 
equations of the Bondi problem. In Section 2 the field equations as well as a brief description of their relevant aspect 
are presented. The Galerkin method is applied to the Bondi problem (Section 3) whose the most important step is the 
choice of convenient basis function for the metric functions. In Section 4 we have presented the following code tests: 
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numerical and exact evolution of weak gravitational waves whose exact solution is described by the linearized Bondi 
equations, verification of consistency relations valid at the future null infinity, and the global energy conservation as 
derived from the Bondi formula of mass loss. Finally, in Section 5 we discuss some perspectives of the present work. 

II. THE BONDI PROBLEM 

The metric for axisymmetric and asymptotically flat spacetimes corresponding to the Bondi problem^ takes the 
form 



ds" = ( -e^" - UW^^ ] du" + 2e^''dudr + 2Ur^e^^dud0 
r 



(1) 



where u is the retarded time coordinate for which the outgoing null cones are denoted by M=constant. The radial 
coordinate r is chosen so that the surfaces of constant {u, r) have area 47rr^; the angular coordinates {9, 0) are constant 
along the outgoing null geodesies. The functions 7, (3, U and V depend on the coordinates w, r, and satisfy the 
vacuum field equations R^^^u = organized in three hypersurface equations and one evolution equation, respectively, 
written as 



(sin S~i) r-e , o 



(2) 
(3) 



Vr = -\r^e'iy-P){U,rf + ^^^^^^ + e'^P-'^) 
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(4) 



(5) 



The subscripts r, u and 6 denote derivatives with respect to these coordinates. Notice that the evolution equation is 
the only one that has a derivative with respect to u, while the hypersurface equations contain only derivatives in the 
null hypersurfaces u =constant. In order to guarantee regularity of the spacetime at the origin r = we have 



7-0(r2), ^_0(^4)^ [/-0(r), V'^r + Oir^). 
Also by imposing smoothness of the axis it is necessary that 



(6) 
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(7) 



be continuous at 6 — 0, tt. 

As pointed out by Bondi et al[l| the above set of equations has a nice hierarchical structure typical of the charac- 
teristic formulation of the Einstein's equations Once the initial data is given, in other words when the function 7 
is specified on the initial null surface u = Uq, the constraint equations ([2|), ([3|) and (j4|) determine the metric functions 
/3, U and V (modulo integration constants) on the null surface u = uq. From these results, the evolution equation 
(|5|) provides "f^u on m = uq, and consequently allows the determination of 7 on the next null surface u = ui, and the 
whole cycle repeats providing the evolution of the spacetime. 
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III. THE GALERKIN APPROACH 



An alternative to the usual strategy adopted in finite difference technique of representing a function by numbers in 
a grid, the Galerkin method or any other spectral method @ assume that the function under consideration can be 
well approximated by a series of a suitable set of basis formed by known analytical functions, as for instance, Legendre 
or Chebyshev polynomials. According to the Galerkin method, the main criterion to choose the basis functions is that 
each component must satisfy the boundary conditions. Thus, let us construct a convenient Galerkin decomposition for 
the metric function 7 using, as our guides, the boundary condition at the origin given by Eq. ([6]), and the regularity 
of 7/ sin^ 6* at = 0, tt. We have come up with the the following decomposition 



j{u,r,x) = (l-x^) J2T. ak,{uM'\r)PAx), (8) 

j=0 k=Q 

where for the sake of convenience we have introduced x = cos 9. Nr and are the truncation orders of the radial 
and angular expansions, respectively; akj{u), k = 0, 1, .., TVj., j = 0, l,..,7Va; are the modal coefficients that depend 
on the retarded time u, and the Legendre polynomials Pj{x) were selected as the angular basis functions due to the 

regularity of 7/ sin^ 6* = 7/(1 — x"^) at x = ±1. Concerning the radial expansion ^'|7''(r) is the basis function of order 
k constructed using a suitable linear combination of the rational Chebyshev polynomials [l^ (see Appendix) in order 
to assure the fulfilment of the boundary condition at the origin, that is ■^'i^\r) ~ ©(r^) near r = for all k. In 

addition the choice of ^'^7 ''(?') provides the correct asymptotic expression for 7 at future null infinity by expanding 
Eq. ([5]) in inverse powers of the radial coordinate, or 



7 = if(ti,x) + ^^ + 0(r-2), (9) 
r 

where the functions K{u,x),c{u^x) are expressed in terms of the modal coefficients akj{u) and the Legendre poly- 
nomials. It is important to emphasize that a global representation of the metric function 7(7/, r, x) as given by Eq. 
([8]), and the other metric functions (see below) allow us to access directly important quantities defined at future null 
infinity that determines the amount of mass and momentum radiated away by gravitational waves. 

Now, the next step is to provide appropriate decompositions for the metric functions f3, U and V . Before establishing 
these decompositions, we follow Ref. [13| and introduce the function S{u,r,x) related to y by ^ = r + r^S. As a 
consequence, the behavior of S near the origin is the same of U and therefore the radial basis for both decompositions 
will be the same. The Galerkin decompositions for the metric functions /3, U and S are 



(3iu,r,x) = il-x'Y c,,{u)^'^>{r)P,{x) 

j-Q k=0 



N!: 



Uiu,r,x) - 6,,(u)vE'i^)(r)P,(a:) 



j=Q k=0 



n; n: 



5(^,r,x) = ^^..,(^)vI/f'(r)P,(x). 

j=0 k=0 



(10) 

(11) 
(12) 



Notice that the angular basis functions are the Legendre polynomials, and each radial basis functions (cf. Appendix) 
must satisfy the boundary conditions given by Eqs. ^ and ([7]) and must provide the correct asymptotic forms for the 
metric. The truncation orders for each Galerkin decomposition are indicated by {N^ , Nx)i {^r ^ ^x) ^'^d {-^r ^ ^x)- 
The corresponding asymptotic expressions for the above metric functions are 



P = H{u,x)+0{r-^) 



(13) 
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U = L{u,x) + 0{r-'^) (14) 

y = So{u,x) H — h ©(r ), (15) 

where it can be shown that the 5*0 — —d{Ly/l — x'^)/dx, S^i is related to the functions H,K,l{^, and M{u,x) 
is known as the Bondi mass aspect [2]. As we are going to see these asymptotic expressions wih be of fundamental 
importance for the evaluation of the Bondi mass and the news function. In the standard Bondi coordinates K — L — 
H = 2]. 

The next step towards the construction of our algorithm to solve numerically the Bondi equations is to substitute 
the decompositions of the metric functions into the field equations. We begin by inserting Eqs. ([8]) and (fTO|) into the 
hypersurface Eq. ([2]) yielding the following residual equation 

. (n^ y 

Res^(«,r,x)=^5]c,,(^^)vI/gp,■(a;)--r ^ E H^S^:' (^) ' (^6) 

j=Q k=0 \j=0 fc=0 J 

Since the decompositions for /3 and 7 are indeed approximations, the residual equation does not vanish exactly. 
Nevertheless, it is expected that the residual approaches to zero as the truncation orders are increased, or Res^ 
as Nr, Nx, , — > 00. According to the basic procedure of the Galerkin method the scalar products, or the 
projections of the residual equation onto all radial and angular basis functions must vanish. These operations are 
schematically indicated by 



X P,{x) f"^ )dx^O, (17) 
^/r{r + 1)/ 

where k = 0,1.., and j — 0,1.., N,^, and the denominator is the weight function of the internal product of 
Chebyshev polynomials. After the integrations we obtain a set of {N!^ + 1) x {N^^ + 1) algebraic equations connecting 
the modal coefficients Ckj with akj . Therefore, the hypersurface Eq. ([2]) is transformed into an algebraic system of 
equations. 

Repeating the same procedure taking into account the hypersurface equations ([3]) and ^ for the functions U and 
S, respectively, we ended up with two set of algebraic equations: the first set connects the modal coefficients bkj with 
the modal coefficients Ukj and Ckj, and the second set that relates the modal coefficients Skj with Ukj, bkj and c^j. 
Finally, after imposing that the projections of the residual equation corresponding to the evolution equation into all 
'^l?^ and Pj{x), a system of (Nx + 1) x (TV^ + 1) nonlinear ordinary differential equations for the modal coefficients 
Ofej emerges. These equations are schematically represented as 



akj{u) = Fkj{aijn,blm,Clm,Slm), (18) 

where dot means derivative with respect to u and k — 0,1.., Nr, j — 0,1.., N^. Therefore the Bondi problem is 
reduced to a system of ordinary differential equations, or a dynamical system along with three sets of algebraic 
constraint equations. In order to integrate the Bondi equations, or equivalently, the dynamical system (|18p . the initial 
conditions aimiuo), 6;m(i*o), Q,„(uo)i s;m(uo) must be given. To accomplish with such a task we have first to determine 
the Galerkin decomposition of the initial data 7o('", x) = 7(1*0, r, x) that will fix the set of aimiuo). The initial values of 
the modal coefficients fe/m(uo) are determined using the first set of algebraic equations obtained from the projections 
indicated by (|17p : considering these results, the second and third sets of algebraic equations determine all Qm(uo) 
and s/m(uo), respectively. In this way, using these initial modal coefficients the dynamical system (jlSp determines the 
modal coefficients aim in the next step, or in the next hypersurface u + Su, and the whole process repeats, providing, 
in this way the numerical solution of the system. 



IV. CODE TESTS 



We now proceed to present three code tests d, [13, [T3| in order to check the issues of stability, accuracy and 
convergence of our algorithm. The first test consists in to evolve a small amplitude gravitational wave and compare 
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it with the corresponding exact solution known from the hnearized field equations. The second test is to check two 
consistency relations at the spatial null infinity obtained from the field equations. In the last test we are going to 
check the global energy conservation given by the Bondi formula that connects the loss of mass with the total energy 
emitted by gravitational waves, where this last quantity is evaluated using the news function. 

In order to carry out these tests we have derived a simplified version of the dynamical system (jlSp . where the 
nonlinear effects are restricted to be of second order in the following way. We have chosen the metric function 7 small 
in the sense of |7max| ^ 10"'^, and with associated modal coefficients \akj \ ^ 10"'^. From the hypersurface equations 
it follows that 0{f3) ~ 0(7^), 0{U) ~ 0{S) ^ C(7), and consequently similar relations expressed in terms of the 
modal coefficients akj , bkj , Ckj and Skj arise. Then, we have generated a dynamical system in which the terms smaller 
or of order, say 0{/3'^), 0(7^), were neglected. The main motivation for this procedure is to provide a fast and simple 
algorithm whose accuracy and convergence can be tested without much computational effort, and at the same time 
catching the main features of nonlinear effects. 



A. Linearized gravitational waves 

In the linear regime /3 ~ and ~ r, and the weak disturbances of the spacctime are described by the functions 
7 and U. It can be shown[3, ^ that the linearized equations are equivalent to a flat scalar wave equation □$ = 0, 
where r, x) is related to 7 and ?7[(|, and hence from the exact solution of the wave equation the corresponding 
exact expressions for 7 and U are determined. Following this procedure the quadrupole solution denoted by 7[2] is 
expressed as 



7[2](u,r,x) 



6Ao{l-x'^)r^{u + r + l) 



(u + l)4(7i + 2r + l)3 
where Aq is the free amplitude of the linearized equations. 



(19) 
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FIG. 1: Illustration of the agreement between the exact (continuous line) and approximate (boxes) expressions for the 
quadrupole mode (|19|l . The curves from the top to bottom are evaluated at u = 0,0.1,0.3,0.5, 1.0, 1.5 and 2.5, respectively. 
The truncation orders are Nr — 6, Nx = 4 and — 7, = 5 for the Galerkin decompositions of 7 and U, respectively. 

Now we chose a very small value for the amplitude, say Aq — 5.0 x 10^**, in order that the nonlinear code be able to 
reproduce the linear solution. The initial values of the modal coefficients akj {u — 0) are determined from the above 
initial data after establishing the following residual equation 



Residual = 7[2] (0, r, x) 



E 

k,j=0 



(20) 



where 7 = (1 — 3:^)7, and impose that their projections into the radial and angular basis functions vanishes, or 
-fW P ^-^aV _ n ^ _ n 1 AT ,.„A ^ - n ^ ^V^, In this way, a set of {Nr + 1) x (A^^ + 1) equations 



^Residual, '^'^'Pn{x)J = for m = 0, 1, ...,Nr and n = 0, 1, 
is generated that can be solved for the initial values akj{0). 
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FIG. 2: (a) Decay of the exact shear K{u,x). (b) Log-linear plot of the absolute error \Kcxmct — J^'numorl with the curves from 
top to below correspond to the truncation orders Nr — 3, ~ 3, Nr = 4, A'^^ = 4 and A^,- — 6, A^'^ — 4 and N^, respectively. 

In Fig. 1 we illustrate the excellent agreement between the the exact (continuous line) and the approximate (boxes) 
solutions evaluated at several times for iV^ = 6 and — 4, where for the sake of convenience we have introduced 
the variable y related to the radial coordinate through r = y/{^ — y) such that < y < 1 corresponds to the 
< r < oo. The decay of the gravitational wave shear 7^] evaluated at future null infinity (y = 1), or the function 
K{u,x) = K/{1 — x^) (cf. Eq. ^) is shown in Fig. 2(a), where for the quadrupole mode K — K{u). It can be seen 
from Fig. 2(b) the evolution of l-R'cxact — ^numcrl with increasing truncation orders {Nr, N^). As expected the greater 
is the truncation orders the numerical solution becomes more close to the exact one. 

At this point it is useful to comment the role played by the choice of the truncation orders associated to the 
Galerkin decompositions of the metric functions. The truncation orders NriN^ of the radial and angular expansions 
of 7, respectively, are our seeds for which the corresponding truncation orders , , , ... will be related with 
through the constrained equations (H]) - (51). We have found that the stability of the code, as far as the evolution 
of weak gravitational waves are concerned, is achieved if the following relations between the truncation orders are 
satisfied, > + I and iVf > A'^j, + 2 for the angular part, and A^^^ > Nr, > Nr associated to the radial 
decomposition. For the sake of convenience we have considered in all numerical experiments the following radial and 
truncation orders: A^^^ = A^,. + 2, A^^^ = A^^ + 1, A^;^ = A^r - 1, and A^^ = A^^; + 1, N^ = N^ + 2, N^ = N^ - 1. 



B. Consistency relations 

The second code test will be to verify two consistency conditionsflS^ that arise from the field equations (jH) and ^ 
evaluated at future null infinity. These relations are written as 



7.« 



S + [{1 ~ x'^)U]' = 
ic7' - (1 - x^)Uj' + 2xU^ = 0, 



(21) 
(22) 



where prime indicates derivative with respect to x, and U = U j \J\ — x^ . To illustrate the type of relations we are 
dealing with in applying the Galerkin method, wc have considered, for instance the term 7^u evaluated at y = 1 
(r = 00) that is given by 



(7, 



u)y=\ 



-am + 772;aoi + 2; - - ao2 



ao.3 + ••• + 7^10 + . 



where the derivatives of the modal coefficients are determined from the dynamical system psp . 



7o(r,x) = 7[2](0,r,a;) 



(u+l)4(u + 2r + l)3 



(23) 



(24) 
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FIG. 3: Relative error S between exact and approximate initial data function. The truncations orders are A',. = 4, Nx = 4, 
A^r = 5, Nx = 4 and A^r = 7, Nx ~ 5. Notice that even for a modest truncation orders Nr = Nx — 4 the fit is quite satisfactory. 
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FIG. 4: Sequence of the log |afej(0)| for the non-vanishing initial modal coefficients for Ao = 0.025 and Nr — 7, Nx = 5, where 
the horizontal axes denotes the order of appearance of the modal coefficients (for instance, aoo(O) — > l,aoi(0) — > 2, etc). 

where 7[2](0,r, x) is the quadrupole mode given by Eq. (jl9p . This initial data can be understood as the quadrupole 
mode plus a self-interaction term (note that the term inside the parenthesis is the quadrupole mode with Aq = 1). 
As shown in Fig. 3, the Galerkin method reproduces quite well the above initial data as indicated by the plots of the 
relative error between the exact and the approximated initial data, 

7o(r, x)-il- x^) Efe,,=o a.,(0)*i^^(r)P,-(x) 

^= 7\ 

loir, x) 

evaluated using distinct truncation orders. 

We determine the initial modal coefficients a^,j(0) from the initial data (IM)) by choosing Aq = 0.025 (see Fig. 4). 
The numerical integration of the dynamical system (|18p was performed using a fourth-order Runge-Kutta method 
in which the three sets of algebraic equations connecting the modal coefficients were taken into account. In this 
way, the evolution of all modal coefficients was obtained and consequently the consistency relations could be verified 
numerically. Fig. 5(a) depicts the log-linear plots of \S+ [(1 — a;^)[/]'| evaluated at m = 0.005 versus x for the following 
truncation orders: A^^ = 4, = 4, iV^ = 5, — 4 and Nr — 7, — 5. As it can be seen the increase of the 
truncation orders implies that this residual equation tends rapidly to zero. However, one could argue that due to the 
fact all the modal coefficients are sufficiently small, that is of order < 10"'^ (cf. Fig. 4), each term we have combined 
are indeed very small producing the curves of Fig. 5(a). For this reason we have plotted in Fig. 5(b) each term, S 
and [(1 — x'^)Uy evaluated at u = 0.005 for A,. = 5,Nx = 4, where both terms range from identical maximum and 
minimum values about ±0.2, ±0.3 almost canceling each other within the angular domain. 




(25) 
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FIG. 5: (a) Log-linear plots of 15* + [(1 — a;^)f7]'| evaluated at u = 0.005. As indicated we have considered distinct radial and 
angular truncation orders to show the convergence of the numerical code, since the increase of the truncation orders implies 
in IS + [(1 — a;'^)C7]'| 0. (b) The contribution of each term, 5* and [(1 — x'^)U]' for A*',. = 5, = 4 is presented, where they 
almost cancel each other producing the corresponding curve of (a). 



In order to provide a good measure of the error in the numerical approximations of the consistency equations 
([2T|) and ((22|) . we define the L2 norm of the residual equations Residuali — S + [{1 — x'^)tj]' and Residual2 — 

l,u - \U' - (1 - x^)U^' + 2xU^ by 




where k — 1,2. It can be noted from Figs. 6(a) and 6(b) the evolution of the L2 norm associated to each residual 
equation for increasing truncation orders. 




C. Energy conservation 



The last and probably most important test is to check the global energy conservation derived from the Bondi 
formula. Following Gomez et al[6j the expression for global energy conservation in the gauge under consideration is 
given by 
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FIG. 7: (a) Sequence of the initial values of the Bondi mass (cf. Eq. (|3ip l corresponding to the initial data (I24p with Aq = 0.025 
and for (A^^ = 3,iVa; =3), {Nr = 4,iV^ =4), (A^^ = 5,iV^ =4), [Nr = 6,iV^ = 4) and (iV^ = 7,N^ = 5) read from the left 
to the right, respectively. The convergence towards the actual value for the initial mass is indicated, (b) The relative error 
\C{u)\/Mb{0) X 100 in the global energy conservation is evaluated up to m = 1.2, where the dynamical system was integrated 
with stepsize of 10"''. Note that for A*',- = 7,Nx = 5 the energy conservation is attained to about 0.45% accuracy. 



C{u)^Mb{u)~MbM + - 



du 



1 „2ff 



-N^dx = 0, 



(27) 



that is obtained after integrating the Bondi formula. Here Mb{u) is the Bondi, H{u,x) is given by Eq. and N 
is the news function. The conformal uj depends on the gauge we are adopting and arises by connecting the metric of 
the two geometry of a unit sphere in the standard Bondi coordinates with the similar expression in our coordinate 
system, or ds^ — d9% + sin^ 0Bd4>% = uP- ((?^ dff^ + sin^ Oc^'^^ dcj)^) yielding 



2e 



K 



(l + a;)e^ + (l-a;)e- 



where 



-dx. 



The news function N{u,r,x) can be written as pi 



(28) 



(29) 



{cue''")' 



(30) 



The Bondi mass depends directly on the mass aspect M (cf. Eq. (jl5|) ) and the terms arising from the gauge under 
consideration. According to Ref. [l^| the Bondi mass Mb{u) is given by 



Mb{u) = Lo-^ {-M + -e-2^[(l - x^y - Axe' - 2c] - e-^^c' {H' + K'){1 - x^) 
j-i '-2 4 

-e-2^c (i7'2 - 2H'K' - K'^){1 - x^) - ^e-"^c [(1 - x'){H" - K") - 4x(H' + K')] S-dx. 



[(1 - x^)iH" - K") - Ax{H' + K')]y 



(31) 



The standard Bondi frame [2] is characterized by a choice of a coordinate system for which H = K = L = and 
consequently uj = 1. Also, in this frame the Bondi mass depends only on M and Eq. (|27p agrees with the original 
expression for the Bondi formula[2]. As an important remark concerning the Galerkin treatment of the field equations 
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is that the functions K, L and M can be read off directly from the asymptotic expressions for the metric functions 
as presented by Eqs. ®, (fT^ - p^ . 

Before verifying the global energy conservation ((27|) . we first evaluate the initial Bondi mass Mb{0) corresponding to 
the initial data ([M)) with = 0-025 considering several truncation orders. Fig. 7(a) shows a sequence of initial values 
of the Bondi mass Mb{uq = 0) signalizing a clear the convergence towards the actual value for the initial Bondi mass. 
The next step consisted in showing the evolution of the relative error in energy conservation, \C{u)\/Mb{Q) x 100, 
with the truncation orders iV^ = 5, = 4 and Nr = 7, = 5. As indicated in Fig. 7(b), the energy conservation is 
attained to 1.3% and 0.45% accuracies, respectively. 



V. DISCUSSION 



In this paper we have constructed an alternative code for the evolution of the Bondi problem based on the Galerkin 
method. The important feature of a satisfactory accuracy with moderate computational resources, along with very 
fast convergence were indeed observed through the code tests of Section 4. As a matter of fact, the highest truncation 
orders we have tested (Nr — 7, = 5) resulted in a dynamical system of 48 equations, which can be considered 
reasonable in the realm of large numerical simulation using Galerkin method [l2|. It is important to mention that 
the radial and angular truncations associated to each Galerkin expansion of the metric functions /3, U and S (cf. 
Eqs. (jlOp . (jlip and Hl2\i ) can be chosen freely starting from their corresponding minimum values to comply with the 
requirements of stability of weak gravitational waves. In deriving the dynamical system we have set for the sake of 
convenience {N^ = N,. - 1, TVf = TV^, - 1), {N^ = Nr + 2,N^ = + 1) and {N^ = Nr + 1, = + 2). 

Possibly the evaluation of a large number of spatial integrations as demanded by the projections of the field 
equations onto the angular and radial basis functions might be seen as the major computational drawback of the 
present scheme. In particular, the increase of the truncation orders produces a dramatic growth of the number of such 
spatial integrations due to the nonlinear terms present in the field equations. In order to circumvent this potential 
problem we have developed efficient symbolic routines in Maple for these operations, which proved to be low-cost 
computationally. We have also constructed a 4th-order Runge-Kutta integrator to evolve the dynamical system (|18p . 
which is adapted to solve numerically the constraint equations derived from the hypersurface equations at each time 
step of integration. 

There are two important further directions of this work that are currently in progress. The first is to extend the 
present code to encompass the full nonlinearities of the field equations. The second is to consider the field equations 
in the gauge K — L = that can be done after a change in the coordinates u, 6 followed by a rescale of the radial 
coordinate 2] . In this case distinct radial basis functions for 7, U and V have been already determined, and this 
particular gauge will provide a simpler the evaluation of the mass function and the news function. 



APPENDIX A: RADIAL BASIS FUNCTIONS 

We shall present the radial basis functions '^^^\ and which are constructed using linear combinations of 
the rational Chebyshev polynomials [lOj TLk{r) defined in the semi-infinite range < r < 00. First it is necessary to 
define an auxiliary function ipk{i") 

Mr)=TLk+iir)+TLk{r). (Al) 
The radial basis functions for the decompositions of 7, U ot S and /3 are given by 

^rir) = lMr) (A2) 

- \ + ' (A3) 

t(/3)/ N (2fc3 + 14fc2 + 25fc + 9) , 3{k + l) , ,x 

^ -4^'^'^ + 4(7 + 6fc + k-)ik + 3) ^^^^^'^^ - 4(F + 6fc + 7){k + 3)(2fc + 7) ^^+^^^^ " 
(fc + 1)(4A;^ + 38P + iisfc^ + 137fc + 48) (fc + 1)(3 + 2fc)(2 + 4fc -f- P)(fc -f 2) 

4(fc2 + 6fc + 7)(A: + 3)(fc + 4)(2fc + 7) ^^''^ 4(7 + 6fc + P)(fc + 3)(fc + 4)(2A: + 7) ^^""^ ^ 
which reproduces correctly the boundary conditions at r = and r — s- cx). 
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